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Do Newtonian large-scale structure simulations fail to include relativistic effects? 
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The Newtonian simulations describing the formation of large-scale structures do not include rela¬ 
tivistic effects. A new approach to this problem is proposed, which consists of splitting the Hawking- 
Hayward quasi-local energy of a closed spacelike 2-surface into a “Newtonian” part due to local 
perturbations and a “relativistic” part due to the cosmology. It is found that the Newtonian part 
dominates over the relativistic one as time evolves, lending support to the validity of Newtonian 
simulations. 
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I. INTRODUCTION 

The study of primordial density fluctuations, of their 
growth history, and of the large scale structures that they 
seed is one of the focuses of modern cosmology. Density 
fluctuations and the temperature fluctuations that they 
induce in the cosmic microwave background are one of 
the foundational observational tools of theoretical cos¬ 
mology. The growth of large scale structures from their 
initial seeds depends on the theory of gravity, as well as 
the cosmological parameters and there are proposals to 
test the theory of gravity using cosmology [Ij (includ¬ 
ing satellite missions) in addition to studying dark en¬ 
ergy and early physics such as baryon acoustic oscilla¬ 
tions. Several redshift surveys have mapped large scale 
structures in the sky in the last few decades. The ob¬ 
servational results are compared with theoretical predic¬ 
tions, but there is a potentially very serious problem here, 
which has received proper attention only recently M- 
Most predictions for the generation and dynamical evolu¬ 
tion of large scale structures, which are based on massive 
Af-body simulations, are done in the Newtonian limit in 
the dust-dominated era of our universe. While they do, 
of course, account for the expansion of the universe, these 
simulations are essentially Newtonian. A priori, Newto¬ 
nian physics is appropriate on small scales, but not on 
scales comparable to the Hubble radius H~^. The in¬ 
creasing computational power has led to a corresponding 
increase in the number of particles N which now reaches 
10® —10^® and in the size of the box used in these simula¬ 
tions, now spanning 0.5—3 Gpc. Since at redshift z ~ 100 
it is ^ 1.5 Gpc, the size of the box easily exceeds 
H~^ in current simulations 0. One must therefore worry 
about including relativistic effects in these simulations. 
These include special-relativistic (velocity) and general- 
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relativistic (cosmological) effects. Naively, one would like 
to superpose the purely local Newtonian gravitational 
field to a purely cosmological “background” field de¬ 
scribed by unperturbed Friedmann-Lemaitre-Robertson- 
Walker (FLRW) space. This is best done using linear 
perturbation theory, if one wants to restrict to first order 
effects in the local potentials. Special-relativistic effects 
are negligible to first order because the peculiar velocities 
of density perturbations are small M- There is a large 
literature on the effects of the cosmological expansion on 
local systems (see Ref. @ for a recent review), while the 
effects of local systems on cosmology through backreac- 
tion have also been the subject of recent debates [1]. The 
recent study in Q assumes that general relativity cor¬ 
rectly describes gravity, that velocity perturbations can 
be neglected, and that first order post-Newtonian poten¬ 
tials ^/>Ar,(^Ar describe the perturbations to FLRW space. 
We follow the same assumptions here. 

There are two ways of dealing with the perturbations of 
FLRW space: either one fixes a gauge (as done in 0 for 
this specific problem), or one adopts a gauge-invariant 
formalism (as done in i). The two recent studies of 
0 and 0 using these different approaches essentially 
agree on the results. The Newtonian simulations do con¬ 
sider the wrong equations for the local potential, the den¬ 
sity perturbations used in the simulations should be cor¬ 
rected, and the initial displacement of particles in the 
simulations is incorrect 0. However, this does not mat¬ 
ter much: in the gauge-dependent discussion the missing 
terms in the equations cancel out so that the potential 
4’sim. in the simulations is computed correctly even at 
large scales. As a result, the N dark matter particles are 
displaced correctly [3. A physical explanation for this 
occurrence, given in Q], is that the scale of the sound 
horizon {i.e., the length travelled by dark matter parti¬ 
cles since the big bang due to their peculiar velocities) is 
much less than the Hubble radius H~^. However, these 
calculations and interpretation are obtained in the New¬ 
tonian conformal gauge and are gauge-dependent. The 
gauge-invariant treatment of Q computes the correct 
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equations to second order in the potentials and includes 
peculiar motions. Using a scheme developed to analyze 
the backreaction of small scale inhomogeneities on large 
scale dynamics, these authors decompose the metric in 
scalar, vector, and tensor perturbations and, by analyz¬ 
ing order by order the Einstein equations for the gauge- 
invariant quantities, build a dictionary between Newto¬ 
nian and relativistic physics. The dictionary is consider¬ 
ably simplified on small scales I <C H~^. The first order 
results agree with those of Ref. Q. As usual, the dis¬ 
advantage of the gauge-invariant treatment is that it is 
much heavier and physically less intuitive and transpar¬ 
ent than the description obtained by fixing the gauge. 
In this paper we propose an approach which is gauge- 
invariant yet physically intuitive and transparent. There 
are no miracles, though: the physical intuition preserv¬ 
ing gauge-invariance is obtained at the price of simpli¬ 
fying (perhaps over-simplifying) the problem. However 
we believe that the approach is an interesting alternative 
and that, since all our knowledge of large scale structure 
is based on A^-body simulations, the limits of validity 
of Newtonian equations in these simulations should be 
scrutinized very closely with a variety of approaches. 


II. A TOY MODEL 

Preserving both the advantage of gauge-invariance and 
that of physical transparency and intuition, we propose a 
toy model in general relativity. This toy model is clearly 
oversimplified and unrealistic but it serves as a good in¬ 
troduction to our approach to the problem and will be 
generalized to realistic situations in the following section. 
The toy model is based on a single, spherically symmetric 
perturbation. 


A. Preliminaries 


and its gradient is 

VaR = 5 ar (1 - ((iAT - r(j)'j^) 


( 4 ) 


to first order in the potential 4>n- 

In general relativity and in the presence of spheri¬ 
cal symmetry, the physical mass-energy contained in¬ 
side a sphere of symmetry is the Misner-Sharp-Hernandez 
(MSH) mass Mmsh 0 ■ It is defined in terms of the areal 
radius R of the enclosing sphere by 0 

1 - = V'^RVcR, (5) 

R 

which gives, for the line element o, 


2Mmsh 

R 


g‘^^VaRVbR 

(1 - 24>n)~^ {I - 4>n - 

1 - 2r,^)v ■ (6) 


Using Eq. (|3]) and again discarding higher order terms, 
this can be re-expressed as the differential equation 


i?2 


d(i)N 

~dR 


ATmsh ■ 


( 7 ) 


Assuming, for sufficiently large i?, that Mmsh — con¬ 
stant, the solution to m as a power series in R ^ begins 
as 


H 

Therefore, for sufficiently large R we have 

Mmsh w 


4’n = — 


rVl - 24>n 


( 8 ) 

( 9 ) 


where we denote the Misner-Sharp-Hernandez mass sim¬ 
ply with m (this change of notation is completely unnec¬ 
essary here but will avoid confusion in the next section). 


We begin from the post-Newtonian, asymptotically 
flat, approximation of the spacetime metric in isotropic 
coordinates 

ds^ = — (1 -I- 2^Ar) djf + (1 — 2(/)jv) -|- r^dil,^2)^ ) 

( 1 ) 

where dV,^ 2 ) — dO^ + sin^ 9 d(p^ is the metric on the unit 
2-sphere. Assuming a single perturbation with spherical 
symmetry, it is 

i’N = i’N{r), (j)N = 4>N{r) (2) 

(we do not specify the form of the potentials and ipN 
yet). In general relativity based on the Einstein-Hilbert 
action, the two potentials coincide, '(/'w = 4>n- 
The areal radiu^ of this spacetime is given by 

R(?') == r (1 - (()Ar(r))-I-e> ((()^) (3) 


B. The spherically symmetric toy model 


The perturbed FLRW metric that we want to consider 
is 


ds^ 




- (1 -h 2'0Ar) dt]'^ 

-I- (1 - 2(j>N) (dr'^ + 


( 10 ) 


in the conformal Newtonian (or diagonal) gauge, where 
we assume that there is only one spherically symmetric 
perturbation: 


i’N=i’N{r), (j)N=4>N{r). (11) 


^ Recall that the areal radius in a spherically symetric spacetime 


is that function R for which SO(3) symmetric surfaces have area 
equal to 
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Here 77 is the conformal time of the unperturbed FLRW 
space related to the comoving time t hy dt = adr]. 
The metric OT is obtained from the post-Newtonian 
metric o by means of the conformal transformation 
dab —>■ dab = ^^dab with conformal factor H = 0 ( 77 ). 

The idea of our toy model is to introduce an effective 
potential $ 7 v = —Mmsr/R and to split it into two contri¬ 
butions which can be compared, one due to the spherical 
perturbation and the other describing the cosmological 
“background”, as in 

= $ 7 v(local) -|- $ AT (cosmological). ( 12 ) 

A priori such an effective potential seems meaningless 
and the decomposition seems highly gauge-dependent, 
but this is not the case, as we show below. The de¬ 
composition we will demonstrate comes from a geometric 
property of the Misner-Sharp-Hernandez mass. 

Due to this decomposition property, and because the 
metric (HI is conformal to the post-Newtonian metric 
o, we are interested in the transformation property of 
the Misner-Sharp-Hernandez mass Mmsh -^msh un¬ 
der conformal transformations, which is easy to derive 
(see Ref. @) 


Mmsh = HMmsh - — VMlVeH - . (13) 

It follows that the Misner-Sharp-Hernandez mass of the 
spacetime m is given by 


Mmsh = am ■ 


= am + 


~ am + 


{1 - a'H'^r^ 

2(1+27/7^) 

n^R^ 

2{l + 2ijN)a‘^ 
oj2 f>3 

^(1-27^.). 


(14) 


where R = a' ja, a prime denotes differentiation with re¬ 
spect to the conformal time 77 of the unperturbed FLRW 
space, and R = a{Tf)-^l — 2(j)ff r = a{r])R is the areal ra¬ 
dius of the metric gab- In terms of the comoving Hubble 
parameter H = a/a = a!/a} = "H/a, we have to first 
order 


Mmsh =ma-\ - - — (1 - 27/;^) . (15) 

The result in Eq. m exhibits a clear decomposition of 
the Misner-Sharp-Hernandez mass into two parts: we de¬ 
note the first term on the right hand side “local part” and 
the second term “cosmological part”. This decomposi¬ 
tion into cosmological and local parts is gauge-invariant, 
as is the full quantity Mmsh- We will demonstrate this 
invariance in a later section when we embed this result 
in the more general case lacking spherical symmetry. 

We can now define the “cosmological effective poten¬ 


tial ”0 


tt2 o2 

(1 - 


(16) 


whose order of magnitude is the square of the physical 
size R of the region of interest in units of the Hubble 
radius and the “local effective potential” 


= — 


ma 

R 


m 

r 


such that the full potential decomposes as 


Mmsh 

R 


$ = 4 *at + . 


(17) 


(18) 


The relative importance of the cosmological and lo¬ 
cal parts of the Misner-Sharp-Hernandez mass (what one 
could call the “degree of non-Newtonianity” of the sys¬ 
tem), is measured by the ratio 


a = 


$Ar 


(1 - 27/7Af) R 

2 ma 

H^R^ 

2ma 

H^R^ 

RSchw 


(19) 


where i?schw = 2ma is the proper Schwarzschild radius 
of the spherical perturbation. The ratio a measures how 
good the Newtonian approximation is: it is good if a <C 1 
and poor otherwise. 

Weakly bound objects are essentially comoving with 
the cosmological “background”, R ~ ar. In the dust- 
dominated era of the universe, during which the forma¬ 
tion of structures takes place, the scale factor evolves as 
a(t) = and we have for weakly bound objects that 


H'^R^ 


4 3 

9 “° 


( 20 ) 


is constant. During this time, the proper Schwarzschild 
radius increases with the scale factor like . Hence 


2r3 1 

9 771 ^2/3 


( 21 ) 


decreases. 

The conclusion achieved with our toy model is that, 
even if numerical simulations of structure formation start 
on scales at which Newtonian physics is not adequate, the 
situation improves with time. For weakly bound objects, 
local physics comes to dominate over the relativistic cos¬ 
mic expansion and this conclusion is gauge-independent 
but also physically intuitive (although the price to pay 
for this simplicity is oversimplification). 

Note that our toy model applies to a perturbed metric 
and fails if there are large primordial black holes. 


^ On small scales R H <I>Ar reduces to the usual Newtonian 
potential —mIR, where m is the Newtonian mass. 
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III. GENERALIZING THE TOY MODEL 

We now wish to relax the unphysical assumption of 
spherical symmetry. The appropriate definition of quasi¬ 
local energy is given by the Hawking-Hayward construct 
Mhh defined for closed spacelike orientable 2-surfaces by 

( 22 ) 

where TZ is the induced Ricci scalar on the 2-surface S, 9± 
and cr^fj are the expansions and shear tensors of a pair of 
null geodesic congruences (outgoing and ingoing from the 
surface S ), a;“ is the projection onto S of the commutator 
of the null normal vectors to S', /i is the volume 2-form 
on S and A is the area of S Q . 

We wish to calculate Mhh for the metric m with 
(pN and ipN general functions of all the coordinates and 
attempt to decompose it into a ‘local’ and a ‘cosmolog¬ 
ical’ piece. Working within general relativity, so that 
ipN = pN, we will do this in two steps: firstly we will cal¬ 
culate Mhh for the post-Newtonian metric m and then 
make use of a known transformation property of Mhh 
under conformal transformations to extend the result to 
the conformally related metric (US but without the as¬ 
sumption of spherical symmetry. We will also calculate 
Mhh directly from the metric m which will allow us 
a finer decomposition of the local part of Mhh into a 
purely Weyl (vacuum) and a local piece coming from the 
perturbation itself. 


A. Non spherically-symmetric post-Newtonian 
metric 


The metric m is a perturbation gab = g^ab + ^9o-b of 
Minkowski spacetime for which = 0 and = 0 
for any 2-surface. Therefore we must have, for the post- 
Newtonian metric, 

^tb = 0 {(I)n) , (23) 

= O , (24) 


hence the squared contributions to Mhh from these ten¬ 
sors are O That is, to first order we have 

Now, from the fully general contracted Gauss equation 
-f e+9_ - a+a'Pp = h^^h^^^Rabcd , (26) 


where hab is the induced metric on the 2-surface S, we 
compute explicitly for the post-Newtonian metric 


n + 9+9. = + 5h^'^')Rabcd 

= hl^^h\^^Rabcd + O {Pi) 


showing this integrand to be of order O (Pn) where we 
have used the fact that the Riemann tensor is first or¬ 
der in the potentials. The “background” inverse induced 
metric is diagonal, 

sin-2 q (28) 

whence 


n + 9+9_ 



(29) 


Using an algebraic software tool such as Maple it is simple 
to compute this Riemann component. After discarding 
all products of pN with any derivative pN and second 
order terms pl, Pn^, one is left with 


= r2 cos 9sm9^ + 2r^ sin^ ^ 


d9 


dr 


2 ■2nd'^^N 


+r sin 9- 


,d'^p 


•N 


902 dp'^ ■ 

Therefore, the final result for Mhh is given by 


(30) 


_ 1 / A f / 2 cos9dpN , 4:dpN 

™ “ SttV 167ris^ Vr2 sin0 d9 ^ r dr 


2 d^pN 2 d^pN \ 

^r2 902 ^ r2sin2 0 dp^^ ) ' ^ ’ 

As a check, assume that pN = pN{r) only. Then the 
first, third, and fourth terms in the bracket of (1311) vanish 
and the second term can be brought out of the integral 
leaving 


Mhh 


1 4 dpN f 

Stt V 167r r dr j 

\ TT / 8r dr 

n2 dpN 

Mmsh 


(32) 


as required (recall the result in Eq. ([7])), where we have 
made use oi R = r O {(Pn) from ([3|) in the second to 
last line. 

In a subsequent section we will be interested in the 
question of gauge invariance for cosmological perturba¬ 
tions. To this end we note that the Riemann tensor is 
decomposed into Ricci and Weyl parts according to m 


Rabcd 9a[c^d\b 9b[c^d\a ^ 9a[c9d\b 4” Cabcd (33) 

and, accordingly, the Hayward-Hawking mass is decom¬ 
posed into Ricci and Weyl parts 


Mmsh = MRicd -I- Mweyi 
Stt \/ IGtt ^ 


\V^Pn + h‘^‘^h'>‘^Cabcd 




(27) 


(34) 
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where is the Laplacian. This is a clean and sharp 
decomposition, which is gauge-independent (we will show 
this below in Sec. IIII C() . 


B. Non spherically-symmetric post-FLRW metric 


1. Calculation of Mhh via conformal transformation 


Using the result ra and a previous result for the 
transformation of the Haw king -Hayward energy under 
conformal transformations [l 2 |, we can compute the 
Hawking-Hayward energy for the post-FLRW metric 
The general result for the conformal transforma¬ 
tion of Mhh is 


Mhh = 


1 

+ 4 i 


VaVfeH 


2Vor2V 





where the metric hat and measure /r, as well as the covari¬ 
ant derivatives Va are those with respect to the confor¬ 
mally related post-Newtonian metric ©• In the special 
case we are interested in, H = fl{r]) only and hence the 
first term inside the integrand vanishes. The second term 
is computed using 


= 

__ ‘^4‘n ,n o 

(l + 20w)2 

~ -2(j)N,n^,n (36) 


where we truncated to first order in the perturbation in 
the last line. The last term in the integrand is simpler 
where we have 


— 

(1 -|- 2(f>]\[) 

= + 2n%ci)N ■ (37) 

Therefore, the final result is 

Mhh = f^Mnn H— (^J -^ j 


2 n 


(38) 


At late times in the dust-dominated era, ^{rj) = Clorf' 
and the right hand side is saturated by the first term, 
which is constant since the second term decays as 77 “^. 
That is, at late cosmological times 


Again, we can perform a check of the general result (l38l) 
under the assumption that (j)M = 4’n{v) only. Writing 
a := Cl for comparison, our formula (13811 simplifies in 
that case to 

= Mhh + (HcfNA) -b 
a 47r 1 

03-^2 

= Mhh H- 2— ~ 2(^jv) , (40) 

which agrees exactly with the formula m written in 
terms of the areal radius in the post-Newtonian space 
R = R/a. 


2 . Direct calculation of Mhh from the metric ng 


An independent calculation of the result (1551) can be 
performed directly by making use of the formula (1551) . 

The metric m is a perturbation of the spatially flat 
FLRW spacetime for which uja vanishes This im¬ 
plies that uJa = O {(j)N) ior the metric (ITUl) and hence the 
squared contribution uJaUJ°‘ can be neglected in the com¬ 
putation of Mhh- Therefore the left hand side of (1551) 
coincides with Mhh up to first order in The right 
hand side of that equality is easy to compute (using Maple 
for example) for the metric g given in (HOI) . 

We wish to decompose the Mhh into ‘local’ and ‘cos¬ 
mological’ contributions, similarly to what we did in the 
spherically symmetric case. Using again the general de¬ 
composition (1331) it can be shown that to first order 







4 ’N,r] 


(41) 


where is the spatial Laplacian in spherical coordinates 
for the post-Newtonian metric gab- Therefore the full 
Hayward Hawking energy is given by 


Mhh 


4TT n 



2 n ' 

(42) 


Each of the terms in (|45|) have a clear physical inter¬ 
pretation. The first two terms are ‘local’, the first of 
gravitational origin being purely Weyl while the second 
being the contribution from the perturbation itself. In¬ 
deed, using the Poisson equation this second term can be 
written as 


^ 2 . f 

gV (j)pf = — l fip^ 


Mhh — f^Mnn ■ 


(39) 


CIR 


(43) 
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showing it to be simply the integrated density of the 
perturbation itself. The third, fourth and fifth terms in 
(1421) are ‘cosmological’, the fifth being the mass energy 
of the cosmological fluid (responsible for the cosmological 
expansion) contained in the sphere of radius R. By the 
Friedman equation this fifth term can be written 


t?3 q2 q 

— -^ = = fIVolpc 

Z \l oTT 


(44) 


where H = 14 is the Hubble parameter and pcos is 
the density of the (unperturbed) cosmic fluid. 

Due to the conformal invariance of the ( i ) Weyl 


tensor Cabc^, we can relate the pure Weyl term to the 
corresponding term computed with respect to the post- 
Newtonian metric by 


ZacZbd/^ _ -Lac-Lbd/^ 

ri fl (-^abcd — ^2 ^ ^ ^abcd 


which leads us to re-derive the result 


R 


Mur = DMhh — ^N,rt + 


14 


14 


(45) 



2 14 ■ 


(46) 


where we recall the result 

This alternative calculation through the use of the 
Weyll tensor allows us the additional refinement of the 
local contribution, which is the only contribution to Mhh 
which survives at late times during dust domination, into 
purely Weyl (vacuum) and purely Ricci (matter) parts as 


Mrr^^ J^il(h--h^^Cabcd) 
= 14Mhh 



These are the first two terms in the full result (H^ . 

This result corroborates the spherically symmetric case 
where we saw that the cosmological part of the quasi¬ 
local energy became negligible at late times in the dust- 
dominated era. 


C. Gauge invariance 

One of the advantages of this work over other ap¬ 
proaches in the literature (for example the works of 
Refs. [l-Q) on the question of the Newtonian approxima¬ 
tion in large scale structure formation is that our results 
are explicitly gauge-invariant, as well as being physically 
transparent. Below we illustrate this gauge invariance for 
both the general and spherically symmetric decomposi¬ 
tions of the quasi-local energy to first order in perturba¬ 
tion theory. 

Firstly, it is known that Mhh, and hence the special 
case Mmsh, is a gauge-invariant quantity. As we have 


already noted, to first order in perturbations the quasi¬ 
local energy Mhh is proportional to the surface integral 
of the scalar quantity 

n + e+e_ - = h‘^^h’^‘^Rabcd. (48) 

This integrand splits naturally into two parts based on 
the gauge-invariant splitting of the Riemann tensor into 
a trace (‘Ricci’) and tracefree (Weyl) component as we 
have already noted: 

R 

Rabcd — 3a[cRd\h 9h[c^d]a 9a[cQd\h 4” ^ahcd ■ (49) 

At late times Mhh converges to precisely the Weyl contri¬ 
bution and an extra piece, which sums the energy density 
of the perturbation itself (in a ‘quasi-local’ fashion, see 
Eq. (H71) l. The Weyl contribution is considered ‘local’ 
since it does not contain information about the cosmic 
expansion (apart from a pre-factor) while the additional 
energy arising from the perturbation, while being part of 
the ‘Ricci’ contribution to Mhh, is clearly local in charac¬ 
ter. The (first order) gauge invariance arises from the fact 
that Eq. (|T71) relates a gauge-invariant quantity, Mhh, 
with two objects, one of which (the purely Weyl contri¬ 
bution) is gauge-invariant. Therefore, to this level of per¬ 
turbation theory, also the second term is gauge-invariant 
and hence the decomposition into local and cosmologi¬ 
cal parts is gauge-invariant at late times. Our conclusion 
that the local part saturates the quasi-local energy at late 
times in a dust-dominated era is also a gauge-invariant 
statement. 


IV. CONCLUSIONS 

Since all our understanding of the formation of large 
scale structures is based on Newtonian simulations, it 
is important to check whether using what are basically 
the wrong equations leads to errors and a misleading 
picture. Confidence based on physical intuition while 
running the simulations is not a substitute for careful 
checking. It is only recently that the works Q and Q 
have taken on this problem and reassured us about the 
validity of the final results of these Newtonian calcula¬ 
tions even if the wrong premises are assumed. In view 
of the importance of this subject, however, it seems that 
more approaches with different techniques should be pur¬ 
sued. The new approach that we present here consists 
of splitting the Hawking-Hayward quasi-local energy of 
a closed, orientable, spacelike 2-surface into a “Newto¬ 
nian” part due entirely to the perturbations of a FLRW 
space (which later become the large scale structures), and 
a “cosmological” part which contains the relativistic ef¬ 
fects. The idea is that, if Newtonian physics is really 
coming to dominate and the relativistic effects remain 
small, then the physical mass-energy of a region of space- 
time should be dominated by a Newtonian contribution 
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while the relativistic effects should give a negligible con¬ 
tribution. We have found and presented the relevant de¬ 
composition of the quasi-local energy; this decomposition 
is gauge-invariant. We have first analyzed a simple toy 
model to illustrate the idea and provide a sense of how 
the calculation proceeds, and then we have moved on to 
realistic situations. For general perturbations, it is not 
trivial to prove that the sought for decomposition of the 
quasi-local energy into Newtonian and relativistic parts 
exists; it follows from the geometric and gauge-invariant 
splitting of the Riemann tensor of the perturbed FLRW 
space into a Ricci part and a Weyl part. At late times 
in the dust era, the Newtonian part of the quasi-local 
mass dominates over the contributions containing the rel¬ 


ativistic effects. The advantage of our approach (limited 
to linear order in the metric perturbations) is that it is 
covariant and gauge-invariant, yet relatively simple to 
follow. Our result provides confidence in the Newtonian 
simulations and agrees with the spirit of the analyses of 
ii, although the methods are quite different. 
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